!************************************************************************
! Gauss-Legendre integration points and weights for the reference triangle
! and tetrahedron. Orders 1, 3, 5 available
! Copyright (c) 2018 Johannes and University of Helsinki
! All rights reserved.
! The MIT License is applied to this software, see LICENSE
module integration_points
   use common
   implicit none

contains

!****************************************************************************80

   subroutine gaussint(P, w, order)
      double precision, dimension(:), allocatable :: P
      double precision, dimension(:), allocatable :: w
      integer, intent(in) :: order

      allocate (P(order), w(order))

      if (order == 1) then
         P = 0.5
         w = 1
      end if

      if (order == 2) then

         P = [0.211324865405187, &
              0.788675134594813]

         w = [0.500000000000000, &
              0.500000000000000]
      end if

      if (order == 5) then

         P = [4.691007703066807d-02, &
              2.307653449471584d-01, &
              5.000000000000000d-01, &
              7.692346550528415d-01, &
              9.530899229693319d-01]

         w = [1.184634425280944d-01, &
              2.393143352496832d-01, &
              2.844444444444447d-01, &
              2.393143352496835d-01, &
              1.184634425280944d-01]

      end if

      if (order == 10) then
         P = [1.304673574141424d-02, &
              6.746831665550779d-02, &
              1.602952158504877d-01, &
              2.833023029353764d-01, &
              4.255628305091844d-01, &
              5.744371694908157d-01, &
              7.166976970646235d-01, &
              8.397047841495122d-01, &
              9.325316833444921d-01, &
              9.869532642585859d-01]

         w = [3.333567215434415d-02, &
              7.472567457529032d-02, &
              1.095431812579910d-01, &
              1.346333596549981d-01, &
              1.477621123573764d-01, &
              1.477621123573765d-01, &
              1.346333596549982d-01, &
              1.095431812579910d-01, &
              7.472567457529045d-02, &
              3.333567215434419d-02]

      end if

      if (order == 20) then
         P = [0.003435700407453, &
              0.018014036361043, &
              0.043882785874337, &
              0.080441514088891, &
              0.126834046769925, &
              0.181973159636742, &
              0.244566499024586, &
              0.313146955642290, &
              0.386107074429177, &
              0.461736739433251, &
              0.538263260566749, &
              0.613892925570822, &
              0.686853044357710, &
              0.755433500975413, &
              0.818026840363258, &
              0.873165953230075, &
              0.919558485911109, &
              0.956117214125663, &
              0.981985963638957, &
              0.996564299592547]

         W = [0.008807003569576, &
              0.020300714900194, &
              0.031336024167054, &
              0.041638370788352, &
              0.050965059908620, &
              0.059097265980759, &
              0.065844319224588, &
              0.071048054659191, &
              0.074586493236302, &
              0.076376693565363, &
              0.076376693565363, &
              0.074586493236301, &
              0.071048054659191, &
              0.065844319224588, &
              0.059097265980759, &
              0.050965059908621, &
              0.041638370788352, &
              0.031336024167054, &
              0.020300714900194, &
              0.008807003569576]

      end if

   end subroutine gaussint

!****************************************************************************80

   subroutine inttri(P, w, order)
      implicit none
      double precision, dimension(:, :), allocatable :: P
      double precision, dimension(:), allocatable :: w
      integer, intent(in) :: order

      integer, dimension(12) :: ind

      ind = (/1, 0, 4, 0, 7, 0, 0, 0, 0, 0, 0, 33/)

      allocate (P(2, ind(order)))
      allocate (W(ind(order)))

      if (order == 1) then
         P(1, 1) = 3.333333333333330d-01
         P(2, 1) = 3.333333333333330d-01

         W(1) = 0.5
      end if

      if (order == 3) then

         P(1, 1) = 3.333333333333330d-01
         P(2, 1) = 3.333333333333330d-01
         P(1, 2) = 0.6
         P(2, 2) = 0.2
         P(1, 3) = 0.2
         P(2, 3) = 0.2
         P(1, 4) = 0.2
         P(2, 4) = 0.6

         W(1) = -2.812500000000000d-01
         W(2) = 2.604166666666665d-01
         W(3) = 2.604166666666665d-01
         W(4) = 2.604166666666665d-01

      end if

      if (order == 5) then
         P(1, 1) = 3.333333333333330d-01
         P(2, 1) = 3.333333333333330d-01
         P(1, 2) = 5.971587178977000d-02
         P(2, 2) = 4.701420641051150d-01
         P(1, 3) = 4.701420641051150d-01
         P(2, 3) = 4.701420641051150d-01
         P(1, 4) = 4.701420641051150d-01
         P(2, 4) = 5.971587178977000d-02
         P(1, 5) = 7.974269853530870d-01
         P(2, 5) = 1.012865073234560d-01
         P(1, 6) = 1.012865073234560d-01
         P(2, 6) = 1.012865073234560d-01
         P(1, 7) = 1.012865073234560d-01
         P(2, 7) = 7.974269853530870d-01

         W(1) = 1.125000000000000d-01
         W(2) = 6.619707639425300d-02
         W(3) = 6.619707639425300d-02
         W(4) = 6.619707639425300d-02
         W(5) = 6.296959027241350d-02
         W(6) = 6.296959027241350d-02
         W(7) = 6.296959027241350d-02
      end if

      if (order == 12) then
         P(1, :) = [2.356522045239000e-02, &
                    4.882173897738050e-01, &
                    4.882173897738050e-01, &
                    1.205512154110790e-01, &
                    4.397243922944600e-01, &
                    4.397243922944600e-01, &
                    4.575792299757680e-01, &
                    2.712103850121160e-01, &
                    2.712103850121160e-01, &
                    7.448477089168281e-01, &
                    1.275761455415860e-01, &
                    1.275761455415860e-01, &
                    9.573652990935790e-01, &
                    2.131735045321000e-02, &
                    2.131735045321000e-02, &
                    1.153434945346980e-01, &
                    2.757132696855140e-01, &
                    6.089432357797880e-01, &
                    2.757132696855140e-01, &
                    6.089432357797880e-01, &
                    1.153434945346980e-01, &
                    2.283833222225700e-02, &
                    2.813255809899400e-01, &
                    6.958360867878030e-01, &
                    2.813255809899400e-01, &
                    6.958360867878030e-01, &
                    2.283833222225700e-02, &
                    2.573405054833000e-02, &
                    1.162519159075970e-01, &
                    8.580140335440730e-01, &
                    1.162519159075970e-01, &
                    8.580140335440730e-01, &
                    2.573405054833000e-02]

         P(2, :) = [4.882173897738050e-01, &
                    4.882173897738050e-01, &
                    2.356522045239000e-02, &
                    4.397243922944600e-01, &
                    4.397243922944600e-01, &
                    1.205512154110790e-01, &
                    2.712103850121160e-01, &
                    2.712103850121160e-01, &
                    4.575792299757680e-01, &
                    1.275761455415860e-01, &
                    1.275761455415860e-01, &
                    7.448477089168281e-01, &
                    2.131735045321000e-02, &
                    2.131735045321000e-02, &
                    9.573652990935790e-01, &
                    2.757132696855140e-01, &
                    6.089432357797880e-01, &
                    1.153434945346980e-01, &
                    1.153434945346980e-01, &
                    2.757132696855140e-01, &
                    6.089432357797880e-01, &
                    2.813255809899400e-01, &
                    6.958360867878030e-01, &
                    2.283833222225700e-02, &
                    2.283833222225700e-02, &
                    2.813255809899400e-01, &
                    6.958360867878030e-01, &
                    1.162519159075970e-01, &
                    8.580140335440730e-01, &
                    2.573405054833000e-02, &
                    2.573405054833000e-02, &
                    1.162519159075970e-01, &
                    8.580140335440730e-01]

         W(:) = [1.286553322022750e-02, &
                 1.286553322022750e-02, &
                 1.286553322022750e-02, &
                 2.184627226901900e-02, &
                 2.184627226901900e-02, &
                 2.184627226901900e-02, &
                 3.142911210894250e-02, &
                 3.142911210894250e-02, &
                 3.142911210894250e-02, &
                 1.739805646535450e-02, &
                 1.739805646535450e-02, &
                 1.739805646535450e-02, &
                 3.083130525779500e-03, &
                 3.083130525779500e-03, &
                 3.083130525779500e-03, &
                 2.018577888319050e-02, &
                 2.018577888319050e-02, &
                 2.018577888319050e-02, &
                 2.018577888319050e-02, &
                 2.018577888319050e-02, &
                 2.018577888319050e-02, &
                 1.117838660115150e-02, &
                 1.117838660115150e-02, &
                 1.117838660115150e-02, &
                 1.117838660115150e-02, &
                 1.117838660115150e-02, &
                 1.117838660115150e-02, &
                 8.658115554329500e-03, &
                 8.658115554329500e-03, &
                 8.658115554329500e-03, &
                 8.658115554329500e-03, &
                 8.658115554329500e-03, &
                 8.658115554329500e-03]

      end if

   end subroutine inttri

!****************************************************************************80

   subroutine inttetra(P, w, order)
      integer, intent(in) :: order
      double precision, dimension(:, :), allocatable :: P
      double precision, dimension(:), allocatable :: W
      integer, dimension(5) :: ind

      ind = (/1, 0, 5, 0, 24/)

      allocate (P(3, ind(order)))
      allocate (W(ind(order)))

      if (order == 1) then

         P(1, 1) = 0.25
         P(2, 1) = 0.25
         P(3, 1) = 0.25

         W(1) = 1.666666666666667d-01

      end if

      if (order == 3) then

         P(1, 1) = 0.25
         P(2, 1) = 0.25
         P(3, 1) = 0.25
         P(1, 2) = 0.5
         P(2, 2) = 1.666666666666667d-01
         P(3, 2) = 1.666666666666667d-01
         P(1, 3) = 1.666666666666667d-01
         P(2, 3) = 1.666666666666667d-01
         P(3, 3) = 1.666666666666667d-01
         P(1, 4) = 1.666666666666667d-01
         P(2, 4) = 1.666666666666667d-01
         P(3, 4) = 0.5
         P(1, 5) = 1.666666666666667d-01
         P(2, 5) = 0.5
         P(3, 5) = 1.666666666666667d-01

         W(1) = -1.333333333333333d-01
         W(2) = 7.499999999999998d-02
         W(3) = 7.499999999999998d-02
         W(4) = 7.499999999999998d-02
         W(5) = 7.499999999999998d-02

      end if

      if (order == 5) then
         P(1, 1) = 3.561913862225449d-01
         P(2, 1) = 2.146028712591517d-01
         P(3, 1) = 2.146028712591517d-01

         P(1, 2) = 2.146028712591517d-01
         P(2, 2) = 2.146028712591517d-01
         P(3, 2) = 2.146028712591517d-01

         P(1, 3) = 2.146028712591517d-01
         P(2, 3) = 2.146028712591517d-01
         P(3, 3) = 3.561913862225449d-01

         P(1, 4) = 2.146028712591517d-01
         P(2, 4) = 3.561913862225449d-01
         P(3, 4) = 2.146028712591517d-01

         P(1, 5) = 8.779781243961660d-01
         P(2, 5) = 4.067395853461134d-02
         P(3, 5) = 4.067395853461134d-02

         P(1, 6) = 4.067395853461134d-02
         P(2, 6) = 4.067395853461134d-02
         P(3, 6) = 4.067395853461134d-02

         P(1, 7) = 4.067395853461134d-02
         P(2, 7) = 4.067395853461134d-02
         P(3, 7) = 8.779781243961660d-01

         P(1, 8) = 4.067395853461134d-02
         P(2, 8) = 8.779781243961660d-01
         P(3, 8) = 4.067395853461134d-02

         P(1, 9) = 3.298632957317306d-02
         P(2, 9) = 3.223378901422757d-01
         P(3, 9) = 3.223378901422757d-01

         P(1, 10) = 3.223378901422757d-01
         P(2, 10) = 3.223378901422757d-01
         P(3, 10) = 3.223378901422757d-01

         P(1, 11) = 3.223378901422757d-01
         P(2, 11) = 3.223378901422757d-01
         P(3, 11) = 3.298632957317306d-02

         P(1, 12) = 3.223378901422757d-01
         P(2, 12) = 3.298632957317306d-02
         P(3, 12) = 3.223378901422757d-01

         P(1, 13) = 2.696723314583159d-01
         P(2, 13) = 6.366100187501753d-02
         P(3, 13) = 6.366100187501753d-02

         P(1, 14) = 6.366100187501753d-02
         P(2, 14) = 2.696723314583159d-01
         P(3, 14) = 6.366100187501753d-02

         P(1, 15) = 6.366100187501753d-02
         P(2, 15) = 6.366100187501753d-02
         P(3, 15) = 2.696723314583159d-01

         P(1, 16) = 6.030056647916491d-01
         P(2, 16) = 6.366100187501753d-02
         P(3, 16) = 6.366100187501753d-02

         P(1, 17) = 6.366100187501753d-02
         P(2, 17) = 6.030056647916491d-01
         P(3, 17) = 6.366100187501753d-02

         P(1, 18) = 6.366100187501753d-02
         P(2, 18) = 6.366100187501753d-02
         P(3, 18) = 6.030056647916491d-01

         P(1, 19) = 6.366100187501753d-02
         P(2, 19) = 2.696723314583159d-01
         P(3, 19) = 6.030056647916491d-01

         P(1, 20) = 2.696723314583159d-01
         P(2, 20) = 6.030056647916491d-01
         P(3, 20) = 6.366100187501753d-02

         P(1, 21) = 6.030056647916491d-01
         P(2, 21) = 6.366100187501753d-02
         P(3, 21) = 2.696723314583159d-01

         P(1, 22) = 6.366100187501753d-02
         P(2, 22) = 6.030056647916491d-01
         P(3, 22) = 2.696723314583159d-01

         P(1, 23) = 2.696723314583159d-01
         P(2, 23) = 6.366100187501753d-02
         P(3, 23) = 6.030056647916491d-01

         P(1, 24) = 6.030056647916491d-01
         P(2, 24) = 2.696723314583159d-01
         P(3, 24) = 6.366100187501753d-02

         W(1) = 6.653791709694644d-03
         W(2) = 6.653791709694644d-03
         W(3) = 6.653791709694644d-03
         W(4) = 6.653791709694644d-03
         W(5) = 1.679535175886776d-03
         W(6) = 1.679535175886776d-03
         W(7) = 1.679535175886776d-03
         W(8) = 1.679535175886776d-03
         W(9) = 9.226196923942397d-03
         W(10) = 9.226196923942397d-03
         W(11) = 9.226196923942397d-03
         W(12) = 9.226196923942397d-03
         W(13) = 8.035714285714281d-03
         W(14) = 8.035714285714281d-03
         W(15) = 8.035714285714281d-03
         W(16) = 8.035714285714281d-03
         W(17) = 8.035714285714281d-03
         W(18) = 8.035714285714281d-03
         W(19) = 8.035714285714281d-03
         W(20) = 8.035714285714281d-03
         W(21) = 8.035714285714281d-03
         W(22) = 8.035714285714281d-03
         W(23) = 8.035714285714281d-03
         W(24) = 8.035714285714281d-03
      end if
   end subroutine inttetra
   
!****************************************************************************80

   subroutine sample_points(P, w, M, N)
      double precision, dimension(:, :), allocatable :: P
      double precision, dimension(:), allocatable :: w
      integer, intent(in) :: M, N

      double precision, dimension(:), allocatable :: P_theta, w_theta
      integer :: i1, i2, j1
      double precision :: phi
      call gaussint(P_theta, w_theta, M)

      allocate (P(3, M*N), w(M*N))

      P_theta = pi*P_theta
      w_theta = pi*w_theta

      j1 = 1
      do i1 = 1, N
         phi = dble(i1 - 1)*2*pi/dble(N)
         do i2 = 1, M
            P(1, j1) = cos(phi)*sin(P_theta(i2))
            P(2, j1) = sin(phi)*sin(P_theta(i2))
            P(3, j1) = cos(P_theta(i2))
            w(j1) = W_theta(i2)/dble(N)*sin(P_theta(i2))*2*pi

            j1 = j1 + 1
         end do
      end do

   end subroutine sample_points

!****************************************************************************80

   subroutine sample_points_cube(P)
      double precision, dimension(:, :), allocatable :: P

      allocate (P(3, 98))

      P(1, :) = [-1.000000000000000, &
                 -0.50000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000]

      P(2, :) = [1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 0.500000000000000, &
                 0.0, &
                 -0.500000000000000, &
                 0.500000000000000, &
                 0.0, &
                 -0.500000000000000, &
                 0.500000000000000, &
                 0.0, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.500000000000000]

      P(3, :) = [-1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 -1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 1.000000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.0, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 -0.500000000000000, &
                 0.0, &
                 0.0, &
                 0.0, &
                 0.500000000000000, &
                 0.500000000000000, &
                 0.500000000000000]

   end subroutine sample_points_cube

end module integration_points
